Ground— state properties of trapped Bose— Fermi 
mixtures: role of exchange— correlation 
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We introduce density-functional-theory for inhomogeneous Bose-Fermi mixtures, derive the as- 
sociated Kohn-Sham equations, and determine the exchange-correlation energy in local density 
approximation. We solve numerically the Kohn-Sham system and determine the boson and fermion 
density distributions and the ground-state energy of a trapped, dilute mixture beyond mean-field 
approximation. The importance of the corrections due to exchange-correlation is discussed by com- 
parison with current experiments; in particular, we investigate the effect of of the repulsive potential 
energy contribution due to exchange-correlation on the stability of the mixture against collapse. 
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I. INTRODUCTION 

The achievement of Bose-Einstein condensation in 
trapped, dilute alkali gases 0| has stimulated a rapidly 
growing activity in the field of ultracold, degenerate 
quantum gases, aimed at a better understanding of fun- 
damental aspects of quantum theory. In particular, re- 
cent experimental progresses have opened the way to the 
fascinating prospect of realizing a BCS transition to su- 
perfluidity in ultracold, trapped Fermi gases. 

Magnetically trapped fermions interact very weakly, as 
their spins are polarized in the direction of the trapping 
magnetic field, so that fermion-fermion s-wave scatter- 
ing is prevented by the Pauli principle. Cooling of the 
fermions to quantum degeneracy can then be efficiently 
achieved by mixing them with ultracold bosons. After 
the process of sympathetic cooling, the final phase of 
the system is a quantum degenerate Bose-Fermi mixture. 
Indeed, such a system has been recently realized experi- 
mentally 0,111113. 

On the theoretical side, dilute Bosc-Fcrmi mixtures 
have been studied both in homogeneous and confined ge- 
ometries. For homogeneous systems, recent work has ad- 
dressed the problem of stability and phase separation Q ; 
the effect of boson-fermion interactions on the dynam- 
ics 0; and the BCS transition induced on the fermions 
by the boson-fermion interactions 0]. The first correc- 
tion to the ground-state energy beyond mean-field ap- 
proximation has been determined analytically for homo- 
geneous systems • This exchange-correlation term can 
be used for trapped systems in local density approxima- 
tion, i.e. when the interaction length scales are much 
smaller than the characteristic sizes of the trapping po- 
tentials. This condition is naturally met in the current 
experiments. Recent numerical work [Tol | confirms the 
analytical findings in the corresponding regime for ho- 
mogeneous systems. 

For trapped systems the theory has been developed 
in mean-field approximation to determine the boson and 
fermion density profiles at zero temperature [Tlj , and the 



related properties of stability against phase separation 
and collapse A mean field approach has been also 
employed to calculate the critical temperature of Bose- 
Einstein condensation in a trapped mixture [T^ . How- 
ever, a description beyond mean field is needed either 
when the interaction parameters are large, or to gain a 
very precise knowledge of the density profiles and the 
related properties of stability. In the present work we 
determine the ground-state energy and the boson and 
fermion density profiles to second order in the boson- 
fermion scattering length for harmonically trapped Bose- 
Fermi mixtures at zero temperature, and determine the 
modification, due to the resulting exchange-correlation 
energy, of the mean-field predictions. 

The plan of the paper is the following. In Section Hll 
we briefiy show how to apply Density Functional The- 
ory (DFT) to inhomogeneous boson-fermion sys- 
tems, and we determine the exchange-correlation energy 
functional via local density approximation (LDA) on the 
ground-state energy functional of homogeneous mixtures 
beyond mean field obtained in Rcf. 0. In Section Uni we 
provide the numerical solution of the coupled, nonlinear 
Kohn-Sham equations for the boson and fermion den- 
sity distributions, and we determine the importance of 
the corrections due to exchange-correlation by compar- 
ing our results with current experiments. In Section llVI 
we discuss the effect of the exchange-correlation energy 
term on the phase diagram of the mixture, especially re- 
garding the onset of collapse for mixtures with attractive 
boson-fermion interaction. 



II. THEORY 

We begin by considering a inhomogeneous, dilute sys- 
tem of interacting bosons and spin-polarized fermions 
with two-body interactions in s-wave scattering ap- 
proximation, so that the interparticle potentials are 
UbbHv - r'l) = gBsSir - r'), UffHt - r'|) = 0, and 
t^BF(|r — r'l) = 9BFS{r — r'). The boson-boson coupling 
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is 5-B-B = '^T^fi^cLBB/nT-B, where ubb is the boson-boson 
s-wavc scattering length, and is the boson mass. 
The boson-fermion couphng reads Qbf ~ 2nh?aBF /itt-r, 
where qb f is the boson-fermion s-wave scattering length, 
and m/j = mBinp / [niB + mp) is the reduced mass {mp 
is the fermion mass). The full Hamiltonian reads: 



H = Tb + Tp + Vb + Vp + Wbb + Wbf: 



(1) 



where 



Tb = - / dr$^(r)^— -$(r); Vb ^ j dr$t(r)T/B$(r), 

fp^- dr¥{v)- ^-fr); Vp ^ / dV^UrWp-^ ir) , 

J 2mp J 

Wbb = \ J y(irdr'$t(r)$t(i.')[/5^$(r')$(r), 

Wbf ^ j j dvdr'^\v)¥{v')UBF-^{Y')^{r). (2) 

Here Tg and Tp denote the boson and fermion kinetic 
energies, Vb{y) and Vp{t) the boson and fermion trap- 
ping potentials, and $(r) and 'l'(r) the boson and fermion 
field operators. 

Let the ground state of the system be \g), and intro- 
duce the ground-state energy =^ {g\H\g), and the 
boson and fermion densities ^^(r) ((7|$'l'(^r)$(r)|g), 

(g|^'l'(r)\I/(r)|g). The Hohenberg-Kohn theo- 
guarantees that, given the interaction poten- 
tials, the ground-state energy depends only on the den- 
sities, i.e. is a functional = EQlnB^np]. The theorem 
was proved originally for Fermi systems, but its gener- 
alization to Bose systems and to Bose-Fermi mixtures is 
straightforward. Determination of the density distribu- 
tions follows by imposing the stationarity conditions: 



nF(rr= 
rem |lj| f 



6Eo[nB,np] \ 



SEo\nB,np] 



Mb 



(3) 



(5nB(r) ' Snp{r) 

where /xb and j.ip arc the boson and fermion chemical 
potentials. 

In general, the functional EQ[nB,np] cannot be deter- 
mined exactly, but we can follow the Kohn-Sham pro- 
cedure to introduce accurate approximations. The 
idea is to map the interacting systems of interest to 
a non-interacting reference system with the same den- 
sity distributions: nB{r) i—> n^^(r) = ^^(r); np{r) i— > 

np {r) = np{r). Uniqueness of the mapping follows from 
the Hohenberg-Kohn theorem, and we find: 

Eo = Tf[nB,np]+Tf[nB,np] 



drVBTiB + / drVpnp 



9BB 
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Qbf J druB-np + £^xc [^b , n^p] , 



(4) 



where the first two terms are the kinetic energies of the 
reference system, the next two terms are the trapping 
energies, and the fifth and sixth term are the mean-field 
part of the interaction energy. The last term includes 
all the contributions to the interaction energy beyond 
mean field due to exchange correlations, and defines the 
exchange-correlation energy functional i?xc['T'B, 'T'f]- If 
Ey,c is neglected altogether, one simply recovers the equa- 
tions of mean-field theory for trapped Bose-Fermi mix- 
tures lillil. 

We now proceed to carry out the full Kohn-Sham 
scheme to determine the ground-state energy, and the 
boson and fermion density profiles beyond mean field. 
In the Kohn-Sham reference system the kinetic parts of 
the energy functional Tg^[nB,np] for the bosons and 
Tp^[nB,np] for the fermions are defined as: 

Tf[nB,np] = -Nb J dM*{r)^^{v), 

Tf[nB.np] = ^Y.J dhrdr^)^^^), (5) 

where Nb and Np are the total numbers of bosons and 
fermions, and the notation (/'(r), is a shorthand for 

the boson and fermion functional orbitals ^[nB,nF](r) 
and ipi[nBTnp\{T) of the non-interacting reference sys- 
tem. Inserting Eqns. jSJ) into Eqn. (0J and carrying 
out the functional derivatives in Eqns. (O, we obtain 
a system of coupled, effective Schrodinger equations for 
the single-particle states that are the desired Kohn-Sham 
equations for a Bose-Fermi system: 
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SE^c 
Sub 

SE^ 
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'ipi=eiipi, (6) 



with ns(r) = iVs|0(r)|2, npir) = j:^^, |^.(r)|2, where 
the sum in np{r) runs over the Np single-particle states 
ipi with lowest energies e^. We now resort to local den- 
sity approximation (LDA) by approximating E'xc with 
an integral over the exchange-correlation energy density 
E!^°"'-{nB{'ic),np{r)) of a homogeneous system taken at 
the -yet unknown- densities ^^(r) and np(r): 



E^c [nB,np]^ / drE^r inB,np) 



(7) 



With this identification, functional derivatives become 
ordinary partial derivatives: 
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dnp 



(8) 



The homogeneous exchange-correlation energy density 
E!^°"^ has been recently determined to second order in 
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the boson- fcrmion scattering length asF via a T-matrix 
approach analog of the Beliaev expansion for a pure Bose 
system , and its expression reads : 



E: 



hom 



{nB,np) 



f{d)kFnFnB, 



(9) 



where kp = {GTr'^npY^^ is the Fermi wave vector, and 
f{5) is a dimensionless function that depends only on 
the boson and fcrmion masses: 

3 + S , 3{l + S)^{l-S)^^^l + 5 



1 



4(5 



1-S' 



(10) 



with S = {ms — mp) /{rriB + mp). Viverit and Giorgini 
have recently shown [l^ that @ is exact in the limit 
kp^B ^ 1, where = 1/ ^/SirnBCiBB is the boson 
healing length. In order of magnitude, the homoge- 
neous densities are np « Np/£^ and ub ~ Nb/(^, 
where i is the characteristic length of the confining po- 
tential. The condition kp^B S> 1 is then equivalent to 
Np > N^J'^{aBB/(-f''^- On the other hand, LDA is cor- 
rect for large Nb and Np, provided that £ ^ clbBtClbf, 
i.e. that the characteristic lengths of the confining poten- 
tials are much larger than the scattering lengths. In cur- 
rent experiments Np « Nb ~ 10'* and aBFl^~ clbb/^^ 
lO^'^, so that the condition fc^^s ^ 1 is well satisfied. 
Moreover, the boson-boson exchange-correlation energy 
is 256h'^aBBn%\/7mB(^^/15mB (see e.g. 0). This is 
much smaller than the exchange-correlation energy (j^ if 

Np » 5AiaBB/aBF)^^'{aBB/e)^^Ha " S)/f{S))N'J\ 
Since aBn/ciBF = 0.13 for the Paris experiment with 
^Li-^Li a and aBB/dBF = 0.28 for the Florence ex- 
periment with ^''K-^^Rb (these are the only two ex- 
periments where qbf has been measured), this condi- 
tion is satisfied as well. Yet other higher order terms 
are due to direct Fermion-Fermion p-wave scattering. 
These terms are at least of order (kpapp)^ , where app is 
the Fermion-Fermion p-wave scattering length, and thus 
certainly negligible against the term we consider. Alto- 
gether, Eq. @ provides the most relevant contribution 
to the exchange-correlation energy for the current ex- 
perimental situations. For more general situations, Eq. 
provides the most relevant contribution beyond mean 
field any time LDA is satisfied, Np is comparable or 
larger than Nb in order of magnitude, and perturba- 
tion theory holds, i.e. kpaBF/i^ « 1; and a sufficiently 
small Bose gas parameter. 

We now consider the Kohn-Sham system © with 
the exchange-correlation energy ^ for spherically 
symmetric, harmonically trapped systems: Vb{tc) = 
(TOBa;|r^)/2, Vp{r) ~ [mpuj'^pr'^) /2. Due to the spheri- 
cal symmetry we can write: 
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(11) 



where Yi„i{<d,^) are the spherical harmonics, and the 
Kohn-Sham equations (jSJ become: 
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nB{r)kF{r) 



Unlir) = eniUni{r) , (12) 



with j dr v?{r) ~ 1, J dr u'f^i(r) ~ 1, where n denotes 
the number of nodes of the radial functions The 
normalized density distributions nB(i") = 47rr^nB(r) and 
np{r) = ATrr'^np{r) are: 



and 



III. 



fiBir) ^ Nbu^t) , (13) 



npir)^ ^ (2; + l)<^(r). (14) 



SOLUTION OF THE KOHN-SHAM 
EQUATIONS 



The above expressions together with Eqns. H12|l de- 
fine a system of coupled nonlinear differential equations. 
The numerical solution is obtained iteratively. We ini- 
tialize nB(r) and ^^(r) to be the Thomas-Fermi den- 
sity distributions with no boson-fermion coupling. We 
then use these as initial densities for Eqns. p2(l . The 
energy eigenvalues are found by a bi-section algorithm, 
iterating the procedure to the desired degree of accu- 
racy. Knowing the states u and Uni, one must deter- 
mine the wave function Uni with lowest energy £„; us- 
ing the fact that £„; grows with n and /. When all 
the occupied Kohn-Sham states are determined, the out- 
put densities are compared to the initial distributions. 
If they are about the same, a self-consistent solution is 
reached, and the procedure ends. If not, one defines a 
convex combination of the initial and output densities, 
«s(F)(^) = (1 - ^) • "b(fT' + ^ • ^'fT^ ^ith < a. < 1, 
and iterates the procedure until convergence is reached 
with the desired degree of accuracy. If Np is large, the 
procedure is very time consuming and limited by a max- 
imum number of nodes that can be included. One then 
adopts a Thomas-Fermi approximation for the fcrmion 
kinetic energy, whenever A''^ > 1000, and finds a poste- 
riori a very good agreement with the single-particle de- 
scription. 

Comparison of our results with current experiments 
can be carried out for those systems whose boson-fermion 
scattering length has been measured. These are the ^Li- 
^Li mixture realized in the Paris experiment 01, and the 
^"K-^^Rb recently realized in the Florence experiment |^. 
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In the Paris experiment with fermionic ^Li and bosonic 
^Li, the measured scattering lengths are ubb = 5. lap 
and asF = 38.0ao, where oq is the Bohr radius. Taking 
ujb &s the unit of frequency, the exchange-correlation en- 
ergy turns out to be 50hujB, whereas the mean- field 
boson-fermion interaction energy is « TA^BHcub- Thus 
only about 0.67% of the interaction energy is due to ex- 
change correlations, it has the same sign of the mean-field 
energy, and the modification of the mean-field density 
profiles is negligible. 

The situation is very different for the mixture of 
fermionic '"'K and bosonic ^^Rb realized in the Florence 
experiment, due to the large and negative boson-fermion 
scattering length giving rise both to a large attractive 
mean-field boson-fermion interaction potential, and to 
a non negligible exchange-correlation potential. The 
latter, being proportional to the square of the boson- 
fermion scattering length, is always repulsive. For this 
experiment, a typical stable configuration is achieved for 
Np ~ 10"*, Nb = 2 X fO^. The boson-boson scattering 
length is aBB = lOOoo, while the boson-fermion scatter- 
ing length aBF ~ — 400ao is measured with an uncer- 
tainty of about 50%. The mean- field interaction energy 
is ~ — 98f65fiw_B, while the exchange-correlation energy 
is Si 6783hujB- Thus the relative correction in the inter- 
action energy is about 7% of the mean- field result, going 
in opposite direction, and leads to a pronounced effect on 
the density profiles. Both the boson and fermion densi- 
ties spread out and decrease substantially at the center 
of the trap with respect to the mean-field prediction, due 
to the repulsive exchange-correlation potential. This ef- 
fect is shown in Figs. ^ and 13 where we show the boson 
and fermion density distributions with and without ex- 
change correlations, calculated with the parameters fixed 
at the values measured in the Florence experiment. At 
the center of the trap the boson and fermion densities 
are reduced, respectively, to about 85% and 78% of the 
mean- field result. 



IV. STABILITY AND COLLAPSE 

In general there are two kinds of instabilities in a bi- 
nary mixture (we do not consider instabilities due to 
fermion pairing): demixing [Tlj and simultaneous col- 
lapse of both the boson and the fermion component |l6j . 
The first can occur if the interaction between the two 
species is repulsive, and implies by definition a minimal 
overlap of the density distributions. In this case we do 
not expect a significant change of the phase diagram by 
repulsive exchange-correlation interactions, but only for 
a small enhancement of the phase separation. 

In the collapse regime, which can occur if the interac- 
tion between the two species is attractive, the situation is 
radically different, as in this case one has indeed a very 
high overlap of the densities in the center of the trap. 
The exchange-correlation interaction, which is always re- 
pulsive to second order in the boson-fermion scattering 
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FIG. 1: The boson density profile for the Florence experiment. 
Dashed line: without exchange correlations; solid line: with 
exchange correlations. Quantities are dimensionless, rescaled 
in units of ^ = (ft/mscjs)"'^''^. 
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FIG. 2: The fermion density profile for the Florence exper- 
iment. Dashed line: without exchange correlations; solid 
line: with exchange correlations. Quantities are dimension- 
less, rescaled in units of ^ = {h/mBi-OBY^^ ■ 

length, opposes the propensity to collapse due to the at- 
tractive mean-field contribution. If the coupling strength 
between the two components of the mixture is sufficiently 
strong, the exchange-correlation can significantly modify 
the phase diagram. 

In Fig. I^lwe provide the mean-field phase diagram of a 
binary boson-fermion mixture, with the physical param- 
eters of the Florence experiment ^3 ■ The plot shows the 
behavior of the critical number of bosons N'^ , i.e. the 
threshold number for the onset of collapse , as a func- 
tion of the number of fermions Np- Collapse occurs at 
any point of the phase plane above the critical curve, 
while the mixture is stable at all points below it. For 
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FIG. 3: The critical number of bosons for the onset of 
collapse as a function of the number of fermions Nf in mean- 
field approximation. 



FIG. 4: The critical number of bosons Ng' for the onset of 
collapse as a function of the number of fermions Nf including 
exchange-correlation. 



low fcrmion numbers Np < 8 x 10^, the critical number 
of bosons begins to grow so fast that to aU practi- 
cal purposes collapse is inhibited. The inversion regime 
between the number of fermions and the critical num- 
ber of bosons takes place at Nf — Ng ~ 5 x 10''^. For 
a typical number of fermions Nf ~ 2 x 10'' one has a 
critical boson number N^ ~ 7 x 10^. The situation in 
mean-field approximation is to be compared with the pre- 
diction obtained by including exchange-correlation. Fig. 
0]shows the same phase diagram as in Fig. |31but with the 
inclusion of exchange-correlation. We clearly see a sig- 
nificant increase in the critical number of the bosons due 
to exchange-correlation. The inversion regime between 
the number of fermions and the critical number of bosons 
takes place at Nf — ~ 1.2 x 10^, and for a typical 
fermion number Np ~ 2 x 10^ the critical boson number 
N^ ~ 1.5 X 10^, i.e. a much larger number of bosons is 
needed to produce a collapse of the fermion component. 
This behavior was qualitatively expected since the effec- 
tive exchange-correlation potentials are always repulsive 
to second order in the boson-fermion scattering length. 

The quantitative difference between the mean-field 
and the exchange-correlation phase diagrams deserves 
some explanatory comments. First of all, the determi- 
nation of the critical line for simultaneous collapse takes 
place in a regime where the numerics is very sensitive to 
small deviations of the input parameters. Thus, when 
a stable solution is not found, this could be ascribed 
either to the fact that the physical coUape regime was 
reached or to an inappropriate numerical precision. How- 
ever, by increasing the numerical precision, computation 
time rapidly increases as well. On the other hand, if a 
stable numerical solution is found, there can certainly be 
no physical collapse. The critical curves we present arc 
then lower bounds on the critical numbers. We remark 



that the mixture is very sensitive to the exact value of the 
boson-fermion scattering length in the collapse regime. 
Since this value is experimentally known with a large 
uncertainty, it would be crucial to determine it with a 
much greater precision. This could be achieved by tun- 
ing the scattering length in order to fit the experimental 
data on the onset of collapse Moreover, for large 

interaction strengths, such as that in the Florence exper- 
iment, the second-order term in the exchange-correlation 
energy might overestimate the effect of stabilization. In 
fact, in these cases, the attractive third-order term could 
possibly give rise to a non negligible contribution, so that 
the mean-field critical line of Fig. (Sjand the second-order 
critical line of Fig. 0] would provide, respectively, a lower 
and an upper bound. The true phase-diagram would 
therefore lie in between the two. A more detailed analy- 
sis than that provided in the present paper requires how- 
ever analytical expressions of the third-order interaction 
energy in powers of kpCLBF, and this is a formidable task, 
because Feynman diagrams containing all possible combi- 
nations of Boson-Fermion and Boson-Boson interactions 
have to be considered. These effects cannot be simply 
determined by resumming restricted classes of equiva- 
lent diagrams. Finally, to go beyond second-order per- 
turbation theory requires, for consistency, to take into 
account interaction processes beyond s-wave, such as p- 
wave scatterings, thus introducing powers of, e.g., the p- 
wave Boson-Fermion scatteting length, and the descrip- 
tion soon becomes exceedingly complex in the framework 
of perturbation theory. Non perturbative methods like 
Monte Carlo simulations would then be desirable to es- 
tablish more accurate results. 

In conclusion, we have introduced the Kohn-Sham 
scheme of DFT for inhomogeneous Bose-Fermi systems 
to determine the ground-state energy and density pro- 
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files to second-order in the boson-fermion scattering 
length. We have compared the theoretical predictions 
with current experiments, discussed the relavance of dif- 
ferent exchange-correlation terms, and investigeted the 
the importance of the exchange-correlation effects for di- 
lute atomic gases. We have shown that they are sub- 
stantial for systems, like ^°K-®'^Rb, with a large attrac- 
tive boson-fermion interaction, especially in the critical 
regime of collapse onset, by comparing the mean-field 



and the exchange-correlation phase diagrams. The DFT 
method outlined here can be in principle extended to 
include higher-order corrections and finite temperature 
effects. 
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